home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
AOL File Library: 2,801 to 2,900
/
aol-file-protocol-4400-2801-to-2900.zip
/
AOLDLs
/
C++ Files Library
/
Graphic Gems I, II & III (C_C++)
/
Graphics Gems C Code.sea
/
GemsIII
/
urot.c
< prev
next >
Wrap
Text File
|
1992-06-16
|
4KB
|
116 lines
/* urot.c */
/* Generates a uniform random rotation */
/* Ken Shoemake, September 1991 */
#include <stdlib.h>
#include <math.h>
#include "GraphicsGems.h"
/* Define an INT32 value to be a 32 bit signed integer */
typedef int INT32;
typedef struct {float x, y, z, w;} Quat;
enum QuatPart {X, Y, Z, W, QuatLen, V=0};
/* * * * * * Utility for quaternion conversion * * * * * */
/** Qt_ToMatrix
* Construct rotation matrix from quaternion (unit or not).
* Assumes matrix is used to multiply row vector on the right:
* vnew = vold mat. Works correctly for right-handed coordinate system
* and right-handed rotations. For column vectors or for left-handed
* coordinate systems, transpose the matrix.
*/
void Qt_ToMatrix(Quat q, Matrix3 *out)
{
double norm = q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w;
double s = (norm > 0.0) ? 2.0/norm : 0.0;
double xs = q.x*s, ys = q.y*s, zs = q.z*s;
double wx = q.w*xs, wy = q.w*ys, wz = q.w*zs,
xx = q.x*xs, xy = q.x*ys, xz = q.x*zs,
yy = q.y*ys, yz = q.y*zs, zz = q.z*zs;
double (*mat)[3] = out->element;
mat[X][X] = 1.0 - (yy + zz); mat[X][Y] = xy + wz; mat[X][Z] = xz - wy;
mat[Y][X] = xy - wz; mat[Y][Y] = 1.0 - (xx + zz); mat[Y][Z] = yz + wx;
mat[Z][X] = xz + wy; mat[Z][Y] = yz - wx; mat[Z][Z] = 1.0 - (xx + yy);
} /* Qt_ToMatrix */
/* * * * * * How to do it using gaussians * * * * * */
/** Qt_RandomG
* Generate uniform random unit quaternion from random seed.
*/
Quat Qt_RandomG(INT32 *argseed)
{
/* This algorithm generates a gaussian deviate for each coordinate, so
* the total effect is to generate a symmetric 4-D gaussian distribution,
* by separability. Projecting onto the surface of the hypersphere gives
* a uniform distribution.
*/
Quat q;
/* uurand generates doubles uniformly distributed between ╨1 and +1 */
/* This linear congruential generator is inline to exploit signed ints */
register INT32 seed = *argseed;
#define uurand() ((seed = (seed+1)*69069)/2147483648.0)
register double x = uurand(), y = uurand();
register double z = uurand(), w = uurand();
register double s1, s2;
double num1, num2, root1, root2, r;
while ((s1 = x*x+y*y) > 1.0) {x = uurand(); y = uurand();}
while ((s2 = z*z+w*w) > 1.0) {z = uurand(); w = uurand();}
/* Now the point (x,y) is uniformly distributed in the unit disk */
/* So is the point (z,w), independently */
num1 = -2*log(s1); num2 = -2*log(s2);
/* Now x*sqrt(num1/s1) is gaussian distributed, using polar method */
/* Similarly for y, z, and w, and all are independent */
r = num1 + num2; /* Sum of squares of four gaussians */
root1 = sqrt((num1/s1)/r); root2 = sqrt((num2/s2)/r);
/* Normalizing onto unit sphere gives uniform unit quaternion */
q.x = x*root1; q.y = y*root1; q.z = z*root2; q.w = w*root2;
*argseed = seed;
#undef uurand
return (q);
} /* Qt_RandomG */
/** M3_RandomRotG
* Generate uniform random rotation matrix from random seed.
*/
void M3_RandomRotG(INT32 *seed, Matrix3 *m)
{
Qt_ToMatrix(Qt_RandomG(seed), m);
} /* M3_RandomRotG */
/* * * * * * How to do it using subgroup algorithm * * * * * */
/** Qt_Random
* Generate uniform random unit quaternion from uniform deviates.
* Each x[i] should vary between 0 and 1.
*/
Quat Qt_Random(double x[3])
{
/* The subgroup algorithm can be condensed to this efficient form.
* Use rotations around z as a subgroup, with coset representatives
* the rotations pointing the z axis in different directions.
*/
Quat q;
register double r1 = sqrt(1.0 - x[0]), r2 = sqrt(x[0]);
register double t1 = PITIMES2*x[1], t2 = PITIMES2*x[2];
register double c1 = cos(t1), s1 = sin(t1);
register double c2 = cos(t2), s2 = sin(t2);
q.x = s1*r1; q.y = c1*r1; q.z = s2*r2; q.w = c2*r2;
return (q);
} /* Qt_Random */
/** M3_RandomRot
* Generate uniform random rotation matrix from uniform deviates.
*/
void M3_RandomRot(double x[3], Matrix3 *m)
{
Qt_ToMatrix(Qt_Random(x), m);
} /* M3_RandomRot */
/* End of urot.c */